source("utils.R")
theme_set(theme_minimal())Irregular Lattices
Dependencies
Until now we have considered the cells to be represented in a point pattern. However, as cells have a shape and area, this might be an oversimplification in some cases. Alternatively, we can rely on the segmentation of indvidual cells that are available for various datasets. The outline of each cell is represented by a polygon and the collection of all cells can be seen as an irregular lattice. Unlike a regular lattice (the representation of spot based spatial tanscriptomics data) the sample areas in an irregular lattice can have different sizes and are not reugularly distributed over the sample space.
For this representation of the cells we will rely on the SpatialFeatureExperiment package. For preprocessing of the dataset we refer the reader to the vignette of the voyager package.
(sfe <- HeNSCLCData())class: SpatialFeatureExperiment
dim: 980 100290
metadata(0):
assays(1): counts
rownames(980): AATK ABL1 ... NegPrb22 NegPrb23
rowData names(3): means vars cv2
colnames(100290): 1_1 1_2 ... 30_4759 30_4760
colData names(17): Area AspectRatio ... nCounts nGenes
reducedDimNames(0):
mainExpName: NULL
altExpNames(0):
spatialCoords names(2) : CenterX_global_px CenterY_global_px
imgData names(1): sample_id
unit: full_res_image_pixels
Geometries:
colGeometries: centroids (POINT), cellSeg (POLYGON)
Graphs:
sample01: # Empty cells
colData(sfe)$is_empty <- colData(sfe)$nCounts < 1
# Select negative control probes
neg_inds <- str_detect(rownames(sfe), "^NegPrb")
# Number of negative control probes
sum(neg_inds)[1] 20colData(sfe)$prop_neg <- colSums(counts(sfe)[neg_inds,])/colData(sfe)$nCounts
# Remove low quality cells
(sfe <- sfe[,!sfe$is_empty & sfe$prop_neg < 0.1])class: SpatialFeatureExperiment
dim: 980 100095
metadata(0):
assays(1): counts
rownames(980): AATK ABL1 ... NegPrb22 NegPrb23
rowData names(3): means vars cv2
colnames(100095): 1_1 1_2 ... 30_4759 30_4760
colData names(19): Area AspectRatio ... is_empty prop_neg
reducedDimNames(0):
mainExpName: NULL
altExpNames(0):
spatialCoords names(2) : CenterX_global_px CenterY_global_px
imgData names(1): sample_id
unit: full_res_image_pixels
Geometries:
colGeometries: centroids (POINT), cellSeg (POLYGON)
Graphs:
sample01: # Calculate count stats
rowData(sfe)$means <- rowMeans(counts(sfe))
rowData(sfe)$vars <- rowVars(counts(sfe))
rowData(sfe)$is_neg <- neg_inds
# log Counts
sfe <- logNormCounts(sfe)In this vignette we will show the metrics related to two marker genes, KRT17 which is a marker of basal cells and TAGLN which is a marker of smooth muscle cells.
plotSpatialFeature(sfe, c("KRT17"),
colGeometryName = "centroids", ncol = 2, scattermore = TRUE) +
theme_void()
plotSpatialFeature(sfe, c("TAGLN"),
colGeometryName = "centroids", ncol = 2, scattermore = TRUE) +
theme_void()
Irregular lattice and neighbourhood matrix
In general a spatial autocorrelation measure has the follwing elements a function \(f(x_i, x_j)\) of the values at locations \(x_i\) and \(x_j\) and a weights matrix \(w_{ij}\). The function \(f\) relates the values \(x_i\) and \(x_j\) and the weights matrix \(w_{ij}\) relates the locations of \(x_i\) and \(x_j\). If \(x_i\) and \(x_j\) are not neighbours the element of the weights matrix is 0. The general the spatial autocorrelation has the of a double sum over all locations \(\sum_i \sum_j f(x_i,x_j) w_{ij}\).
The analysis techniques rely on the a neighborhood matrix, which is a matrix that indicates which cells are neighbors. In the case of the regular lattice, the calculation of the neighborhood matrix is rather straightforward. In the case of the irregular lattice it can be more complicated as the reconstruction of the cells is often not perfect.
Neighbourhood matrix
One of the challenges when working with irregular lattice data is the construction of a neighbourhood graph. The main question is, what to consider as neighbours, as this will affect downstream analyses. Various methods exist to create neighbours, such as contiguitiy based neighbours, graph-based neighbours, distance based neighbours or higher order neighbours. The segmentation of individual cells in notoriously difficult and often the segmentation is not perfect. Therefore we cannot simply rely on neighbours based on the segementation results and the construction of contiguitiy based neighbours is therefore not ideal, also as it would neglect the influence of more distant neighbours.s
References for the construction of the neighbourhood matrix: * https://r-spatial.github.io/spdep/articles/nb.html * https://r-spatial.org/book/14-Areal.html * https://pachterlab.github.io/voyager/articles/vig4_cosmx.html#data-normalization * Getis, Arthur. 2009. “Spatial Weights Matrices.” Geogr. Anal. 41 (4): 404–10. * Zuur, A. F., Ieno, E. N., Smith, G. M., Saveliev, A. A., Mukharamova, S. S., & Zuur, A. F. (2007). Analysis and modelling of lattice data. Analysing Ecological Data, 321-339.
In the following we will use the k-nearest neighbours approach from voyager which relies on spdep. The approach and the number \(k\) is somewhat arbritary.
Thought: if we assume non perfect segmentation we could basically calculate the same statistics also on point pattern derived datasets. We would also use a k-nearest neighbours approach in a “point pattern” dataset. The centroids of the cells as a proxy for the spatial location.
Is there a good methods to dermine a good neighbourhood matrix?
Point pattern as a “event based” dataset, lattice data an “observation” based dataset.
colGraph(sfe, "knn5") <- findSpatialNeighbors(sfe, method = "knearneigh",
dist_type = "idw", k = 5,
style = "W")Global Methods
Global methods give us an overview over the entire study are. The gobal measures summarize the spatial distribution of the cells in a single value. The most common global measures are the Moran’s I and Geary’s C coefficients. These coefficients are based on the spatial weights matrix and the values of the variable of interest.
Global Moran’s I coefficient
The global Moran’s I coefficient is a measure of spatial autocorrelation. It is defined as
\[I = \frac{n}{\sum_i\sum_j w_{ij}} \frac{\sum_i\sum_j w_{ij}(y_i - \hat{y})(y_j - \hat{y})}{\sum_i (y_i - \hat{y})^2}\]
where \(y_i\) is the value of the variable of interest for location \(i\), \(\hat{y}\) is the mean of \(y_i\) and \(w_{ij}\) is the spatial weight between locations \(i\) and \(j\). This value is close to \(0\) for large \(n\). A value higher than \(\mathbb{E}(I) = -1/(n-1)\) indicates spatial auto-correlation. Negative values indicate negative auto-correlation.
Implementation using voyager
calculateMoransI(sfe,
features = c("KRT17", "TAGLN"),
colGraphName = "knn5",
exprs_values = "logcounts"
)DataFrame with 2 rows and 2 columns
moran K
<numeric> <numeric>
KRT17 0.643987 5.09219
TAGLN 0.262379 6.55766sfe <- runUnivariate(sfe,
features = c("KRT17", "TAGLN"),
colGraphName = "knn5",
exprs_values = "logcounts",
type = "moran.mc",
nsim = 200)
res <- rowData(sfe)[c("KRT17", "TAGLN"),]
#value of the metric
res$moran.mc_statistic_sample01[1] 0.6439874 0.2623794#p-value
res$moran.mc_p.value_sample01[1] 0.004975124 0.004975124Implementation using spdep
#create nearest neighbours weights
weights_neighbourhoods <- colGraph(sfe, "knn5")
spdep::moran.test(x = logcounts(sfe)["KRT17",],
listw = weights_neighbourhoods,
randomisation = TRUE)
Moran I test under randomisation
data: logcounts(sfe)["KRT17", ]
weights: weights_neighbourhoods
Moran I statistic standard deviate = 319.35, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
6.439874e-01 -9.990609e-06 4.066523e-06 spdep::moran.test(x = logcounts(sfe)["TAGLN",],
listw = weights_neighbourhoods,
randomisation = TRUE)
Moran I test under randomisation
data: logcounts(sfe)["TAGLN", ]
weights: weights_neighbourhoods
Moran I statistic standard deviate = 130.12, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
2.623794e-01 -9.990609e-06 4.066463e-06 We can also use the moran.mc function to calculate the Moran’s I coefficient. This function uses a Monte Carlo simulation to calculate the p-value. The function is part of the spdep package.
spdep::moran.mc(x = logcounts(sfe)["KRT17",],
listw = weights_neighbourhoods,
nsim = 200)
Monte-Carlo simulation of Moran I
data: logcounts(sfe)["KRT17", ]
weights: weights_neighbourhoods
number of simulations + 1: 201
statistic = 0.64399, observed rank = 201, p-value = 0.004975
alternative hypothesis: greaterWe can see both genes have a positive Moran’s I coefficient and a highly significant p-value. The expected value is \(\mathbb{E}(I) = -1/(n-1)\) which is for large \(N\) close to 0. Positive and significant values indicate clustering of similar like values. It is important to note that this could be both at the high or low end of the values. Negative values indicate clustering of alternating values, i.e. gives a measure of spatial heterogeneity. It is also important to note that the result is dependent on the neighbourhood matrix. Different neighbourhood matrices will give different results. To compare Moran’s I coefficients between different values we need to use the same neighbourhood matrix. Otherwise, we would have to rely on standardized \(z\)-scores.
Global’s Geary’s coefficient
Geary’s \(C\) is a different measure of global autocorrelation and is very closely related to Moran’s \(I\). Howver, it focuses on spatial dissimilarity. Geary’s \(C\) is defined by
\[ C = \frac{(n-1) \sum_i \sum_j w_{ij}(y_i-y_j)^2}{2\sum_i \sum_j w_{ij}\sum_i(y_i-\bar{y})^2} \]
where \(n\) is the number of locations, \(i,j\) are different locations, \(y_i\) is the value of the variable of interest for location \(i\), \(\bar{y}\) is the mean of \(y_i\) and \(w_{ij}\) is the spatial weight between locations \(i\) and \(j\). The interpretation is inversely to Moran’s \(I\). A value less than \(1\) indicates positive auto-correlation, a value more than \(1\) negative auto-correlation.
Implementation using voyager
sfe <- runUnivariate(sfe,
features = c("KRT17", "TAGLN"),
colGraphName = "knn5", nsim = 200,
type = "geary.mc")
res <- rowData(sfe)[c("KRT17", "TAGLN"),]
#value of the metric
res$geary.mc_statistic_sample01[1] 0.3563555 0.7259383#p-value
res$geary.mc_p.value_sample01[1] 0.004975124 0.004975124Implementation using spdep
#create nearest neighbours weights
weights_neighbourhoods <- colGraph(sfe, "knn5")
spdep::geary.test(x = logcounts(sfe)["KRT17",],
listw = weights_neighbourhoods,
randomisation = TRUE)
Geary C test under randomisation
data: logcounts(sfe)["KRT17", ]
weights: weights_neighbourhoods
Geary C statistic standard deviate = 293.94, p-value < 2.2e-16
alternative hypothesis: Expectation greater than statistic
sample estimates:
Geary C statistic Expectation Variance
3.563555e-01 1.000000e+00 4.794895e-06 spdep::geary.test(x = logcounts(sfe)["TAGLN",],
listw = weights_neighbourhoods,
randomisation = TRUE)
Geary C test under randomisation
data: logcounts(sfe)["TAGLN", ]
weights: weights_neighbourhoods
Geary C statistic standard deviate = 121.89, p-value < 2.2e-16
alternative hypothesis: Expectation greater than statistic
sample estimates:
Geary C statistic Expectation Variance
7.259383e-01 1.000000e+00 5.055704e-06 Again the value of Geary’s \(C\) indicates that the genes are spatially auto-correlated.
Local measures
Unlike global measures, that give as an overview over the entire study are local measures give us information about the statistic at each location, i.e. cell. There exist local analogs of Moran’s I and Geary’s C for which the gobal statistc counterpart can be seen as weighted sum of its respective local statistic. These coefficients are based on the spatial weights matrix and the values of the variable of interest.
Local Moran’s I coefficient
The local Moran’s I coefficient is a measure of spatial autocorrelation on each location of interest. It is defined as
\[I_i = \frac{x_i - \bar{x}}{\sum_{k=1}^n(x_k-\bar{x})^2/(n-1)} \sum_{j=1}^n w_{ij}(x_j - \bar{x})\] where the idex \(i\) refers to the location for which the measure is calculated. The interpratation is analog to the gobal Moran’s I. The value of \(I_i\) is close to \(0\) for large \(n\). A value higher than \(\mathbb{E}(I) = -1/(n-1)\) indicates spatial auto-correlation. Smaller values indicate negative auto-correlation. It is important to note that, like for the global counterpart, the value of local Moran’s I could be a result from both the high or low end of the values. Since we measure a great number of locations at the same time, we need to correct for multiple testing. We use the Benjamini-Hochberg procedure to correct for multiple testing.
Implementation using voyager
sfe <- runUnivariate(sfe,
features = c("KRT17", "TAGLN"),
type = "localmoran")
plotLocalResult(sfe, "localmoran",
features = c("KRT17", "TAGLN"), ncol = 2,
colGeometryName = "centroids")
Implementation using spdep
TODO: change this impelemntation to work with gene instead of cell counts
# TODO: add weights matrix, add which sf feature to extract for plotting
#calculate local Moran's I and then correct for multiple testing using Benjamini-Hochberg correction if you want to plot p-values
localplotCell <- function(sfe, var, fun, plotvar, weights_neighbourhoods,
type){
if (type == 'expression'){
#extract the gene counts
x <- logcounts(sfe)[var,]
}
else if (type == 'feature'){
#extract the cell counts
x <- sfe[[var]]
}
else{
stop('type must be either "gene" or "cell"')
}
loc <- do.call(fun, args = list(x= x, listw = weights_neighbourhoods))
#why so ever, 'localG' has a different return structure than 'localmoran'. Thus, this conditional with different indexing
if(fun %in% c('localG')){
loc <- attr(loc, 'internals')
#extract the effect size
locEffect <- loc[,1]
#extract the p-value and adjust for multiple testing
p.val.adj <- loc[,5] |> p.adjust("BH")
}
else if(fun == 'localC_perm'){
p.val.adj <- attr(loc, 'pseudo-p')[,'Pr(z != E(Ci))']
locEffect <- loc
}
else if (fun == 'LOSH'){
locEffect <- loc[,1]
p.val.adj <- c()
}
else{
#extract the effect size
locEffect <- loc[,1]
#extract the p-value and adjust for multiple testing
p.val.adj <- loc[,5] |> p.adjust("BH")
}
#convert into a plain sf object for plotting
sf <- colGeometries(sfe)$cellSeg
sf$locEffect <- locEffect
sf$p.val.adj <- p.val.adj
return(tm_shape(sf) + tm_fill(col = plotvar))
# p <- ggplot() +
# geom_sf(data = sf, aes(fill = plotvar)) +
# scale_y_continuous(breaks = 34:36)
# return(p)
}p <- localplotCell(sfe, 'KRT17', fun = 'localmoran', plotvar = 'locEffect', weights_neighbourhoods = weights_neighbourhoods, type = 'expression')
q <- localplotCell(sfe, 'TAGLN', fun = 'localmoran', plotvar = 'locEffect', weights_neighbourhoods = weights_neighbourhoods, type = 'expression')
tmap_arrange(p,q)
Local Geary’s C coefficient
Similar to local Moran’s I there is a local Geary’s C coefficient. It is defined as
\[C_i = \sum_{j=1}^n w_{ij}(z_i-z_j)^2\]
The interpretation is analog to the global Geary’s C. A value less than \(1\) indicates positive auto-correlation, a value more than \(1\) negative auto-correlation.
In this example we will not plot the local Gear’s C coefficient for gene expression but for features that are associated with the individual cells, e.g. the number oc counts or the number of genes expressed. We will use the colDataUnivariate function to calculate the local Geary’s C coefficient for these features.
Implementation using voyager
sfe <- runUnivariate(sfe,
features = c("KRT17", "TAGLN"),
type = "localmoran")
plotLocalResult(sfe, "localmoran",
features = c("KRT17", "TAGLN"), ncol = 2,
colGeometryName = "centroids")
Implementation using spdep
The same is possible using spdep
p <- localplotCell(sfe, var = 'KRT17', fun = 'localmoran',
plotvar = 'locEffect', weights_neighbourhoods = weights_neighbourhoods,
type = 'expression')
q <- localplotCell(sfe, var = 'TAGLN', fun = 'localmoran',
plotvar = 'locEffect', weights_neighbourhoods = weights_neighbourhoods,
type = 'expression')
tmap_arrange(p,q)
Local Getis-Ord statistic
The local Getis-Ord \(G_i\) statistic quantifies the weighted concentration of points within a radius \(d\) and in a local region \(i\).
\[ G_i(d) = \frac{\sum_{j=1}^n w_{ij}(d)x_j}{\sum_{j=1}^n x_j}, \text{s.t } j \neq i \]
There is a variant of this statistic, \(G_i^*(d)\) which is the same as \(G_i(d)\) except that the contribution when \(j=i\) is included in the term. Importantly, \(G_i(d)\) is scale-invariant but not location-invariant. That means, the subdivision into the \(n\) subregions matters for the computation of the local statistic.
Implementation using voyager
sfe <- runUnivariate(sfe,
features = c("KRT17", "TAGLN"),
type = "localG")
plotLocalResult(sfe, "localG",
features = c("KRT17", "TAGLN"), ncol = 2,
colGeometryName = "centroids")
Implementation using spdep
p <- localplotCell(sfe, var = 'KRT17', fun = 'localG', plotvar = 'locEffect',
weights_neighbourhoods = weights_neighbourhoods,
type = 'expression')
q <- localplotCell(sfe, var = 'TAGLN', fun = 'localG', plotvar = 'locEffect',
weights_neighbourhoods = weights_neighbourhoods,
type = 'expression')
tmap_arrange(p,q)